********************************
** Read in file, add controls **

use 	"./output/cities_expandedpatents", clear

merge m:1 npi_id using "./output/intermediate/controls_43cities"
drop _merge

merge m:1 extra_id using "./output/intermediate/controls_7cities", update
drop _merge

** npi X month-of-year FEs **
egen	idxm = group(gen_id f_m)

** census region X year FEs **
egen	regyr = group(f_yr region)

*=======================*
*= Treatment variables  *
*=======================*

tab 	days_npis if f_myr==tm(1918m9)

gen		daysnpis = days_npis/30 
gen		lnnpis = ln(days_npis)
gen		longnpis = (days_npis>=90) if !mi(days_npis) // Natural break in data, about 61% have short
gen		ldeaths = ln(mag_maxexcessdeaths)

foreach v of varlist daysnpis longnpis lnnpis ldeaths {
	gen 	post_X_`v' 		= `v'*(f_myr>=tm(1918m9))

	gen 	during_X_`v' 	= `v'*( (f_myr>=tm(1918m9)) & (f_myr<=tm(1919m3)) )
	gen 	after_X_`v' 	= `v'*(f_myr>=tm(1919m4))
	gen 	pre_X_`v' 		= `v'*(f_myr<=tm(1917m8))
	
	lab var post_X_`v'		"Post Pandemic $\times$ NPI Length"
	
	lab var pre_X_`v'		"Before Pandemic $\times$ NPI Length"
	lab var during_X_`v' 	"During Pandemic $\times$ NPI Length"
	lab var after_X_`v' 	"After Pandemic $\times$ NPI Length"
}

*======================*
** Summary Statistics **

** Table A.1
estimates clear
eststo : estpost tabstat prate_wtd_inv pat_wtd_inv pat_single_inv pat_multpats_inv pat_wtd_inv_noassg pat_wtd_inv_wassg if tin(1916m1, 1920m12), c(s) s(mean sd min p25 median p75 max)

local 	D = c(current_date)

esttab 	using "$RES/fortable/summary_stats_outcome.csv", replace cells(mean sd min p25 p50 p75 max)
esttab 	using "$RES/fortable/summary_stats_outcome.tex", replace cells(mean sd min p25 p50 p75 max) bookt f
est	 	clear

estimates clear
eststo : estpost tabstat days_npis lnnpis longnpis if f_myr==tm(1918m9), c(s) s(mean sd min p25 median p75 max)

local 	D = c(current_date)

esttab 	using "$RES/fortable/summary_stats_treatment.csv", replace cells(mean sd min p25 p50 p75 max)
esttab 	using "$RES/fortable/summary_stats_treatment.tex", replace cells(mean sd min p25 p50 p75 max) bookt f
est	 	clear

sum pop_int [aw=pop_int] if f_myr==tm(1910m1)
local cities_percent_us = 100*`r(sum_w)'/92410000
di "50 cities contain `cities_percent_us'% of 1910 US Population"

preserve
	collapse (sum) pat_wtd_inv, by(f_yr)
	keep if f_yr==1910
	sum pat_wtd_inv
	local pats1910 = r(mean)
	di "50 cities contain `r(mean)' 1910 US Patent filings, compare to overall number"
restore

** Table A.2
preserve
	keep if f_myr==tm(1910m4)
	keep city days_npis pop_int npi_id
	gen inmarkel = (!mi(npi_id))
	gen source = "Markel et al.\ (2007)" if inmarkel==1
	replace source = "Influenza Archive 2.0" if inmarkel==0
	drop npi_id inmarkel
	sort city
	tostring days_npis, replace
	replace pop_int = round(pop_int)
	recast long pop_int, force
	tostring pop_int, force format(%9.0fc) gen(pop_int_st)
	drop pop_int
	order city days_npis pop_int_st source, first
	* doesn't work on every platform, if challening, instead export delim and them manually put into tex format
	capture dataout, save("$RES/fortable/50cities.tex") tex replace
restore

*=======================*
*= Balance Tables (Table A.3) =*
*=======================*

sort gen_id f_myr 

gen prate_1910 = (L3.prate_all_inv + L2.prate_all_inv + L1.prate_all_inv + prate_all_inv + ///
					F1.prate_all_inv + F2.prate_all_inv + F3.prate_all_inv + F4.prate_all_inv + ///
					F5.prate_all_inv + F6.prate_all_inv + F7.prate_all_inv + F8.prate_all_inv)/12 if f_myr==tm(1910m4)

local setup1 "(mean if longnpis==1) (mean if longnpis==0) (diff longnpis)"
local setup2 "(diff longnpis)"
local dtxt1A "using $RES/fortable/balance_nofe_"
local dtxt1B "using $RES/fortable/balance_wfe_"
local dtxt2 ".tex if f_myr==tm(1910m4), vce(robust) replace starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001)"
 
balancetable `setup1' prate_1910 	`dtxt1A'1`dtxt2' format(%9.2f) sd(par([ ]))
balancetable `setup1' pop_int  		`dtxt1A'2`dtxt2' format(%9.0f) sd(par([ ]))
balancetable `setup1' days_topeakdeaths mag_maxexcessdeaths `dtxt1A'3`dtxt2' format(%9.1f) sd(par([ ]))
balancetable `setup1' sh1910_inschool sh1910_lit sh1910_germandescent sh1910_germanlang sh1910_age20to60 `dtxt1A'4`dtxt2' format(%9.3f) sd(par([ ]))
balancetable `setup1' age1910_census `dtxt1A'5`dtxt2' format(%9.1f) sd(par([ ]))

balancetable `setup2' prate_1910 	`dtxt1B'1`dtxt2' format(%9.2f) fe(regyr)
balancetable `setup2' pop_int  		`dtxt1B'2`dtxt2' format(%9.0f) fe(regyr)
balancetable `setup2' days_topeakdeaths mag_maxexcessdeaths `dtxt1B'3`dtxt2' format(%9.1f) fe(regyr)
balancetable `setup2' sh1910_inschool sh1910_lit sh1910_germandescent sh1910_germanlang sh1910_age20to60 `dtxt1B'4`dtxt2' format(%9.3f) fe(regyr)
balancetable `setup2' age1910_census `dtxt1B'5`dtxt2' format(%9.1f) fe(regyr)

*=======================*
*= Control Trends =*
*=======================*

gen		T = (f_myr - tm(1916m1))/12

gen		sh1910xT_lit = sh1910_lit * T
gen		sh1910xT_inschool = sh1910_inschool * T

gen		sh1910xT_germd = sh1910_germandescent * T
gen		sh1910xT_germl = sh1910_germanlang * T

gen		sh1910xT_age2060 = sh1910_age20to60 * T
gen		sh1910xT_ageave = age1910_census * T

gen		NPIxT = daysnpis*T

lab var sh1910xT_lit 		"Linear Trend in 1910 Literacy Share"
lab var sh1910xT_inschool 	"Linear Trend in 1910 In School Share"

lab var sh1910xT_germd		"Linear Trend in 1910 Share of German Descent"
lab var sh1910xT_germl		"Linear Trend in 1910 Share with German mother tongue"

lab var sh1910xT_age2060	"Linear Trend in 1910 Share age 20-60"
lab var sh1910xT_ageave		"Linear Trend in average age in 1910"

lab var NPIxT 	"Time Trend $\times$ (Days of NPIs)/30"

*=======================*
*= Robustness Table (Table A.4)  ==*
*=======================*
est clear

local treat1 post_X_longnpis
local treat2 pre_X_longnpis during_X_longnpis after_X_longnpis
local outdet "exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)"


eststo: ppmlhdfe pat_wtd_inv `treat1' if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' if tin(1916m1, 1920m12) & in_Markel==1, `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' if tin(1916m1, 1920m12) & in_Markel==1,`outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' post_X_ldeaths if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' pre_X_ldeaths during_X_ldeaths after_X_ldeaths if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' sh1910xT_lit sh1910xT_inschool if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' sh1910xT_lit sh1910xT_inschool if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' sh1910xT_age2060 sh1910xT_ageave if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' sh1910xT_age2060 sh1910xT_ageave if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' sh1910xT_germd sh1910xT_germl if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat2' sh1910xT_germd sh1910xT_germl if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' NPIxT if tin(1916m1, 1920m12), `outdet'
eststo: ppmlhdfe pat_wtd_inv `treat1' i.gen_id#c.T if tin(1916m1, 1920m12), `outdet'

esttab 	using "$RES/fortable/robust_all.tex", booktabs replace cells(b(star fmt(%9.3f)) se(par)) stats(N, fmt(%9.0g))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001)
est 	clear


*=======================*
*= Table A.6 Multi-Inventor Patent Robustness  ==*
*=======================*

foreach v in longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_wtd_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_first_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_all_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	
	estout 	using "$RES/fortable/robustpats_posteffect_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
			stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
			style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

foreach v in longnpis {
	
	local treat pre_X_`v' during_X_`v' after_X_`v'

	eststo: ppmlhdfe pat_wtd_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_first_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_all_inv 		`treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/robustpats_duringafterpre_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
			stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
			style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

*=======================*
*= Table A.7 Population Robustness  ==*
*=======================*

* Baseline

foreach v in longnpis { //daysnpis lnnpis longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr gen_id) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1920m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1923m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1923m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1910m1, 1926m12), exp(pop_int) a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/poprobust_baseline_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}


* No control for population

foreach v in longnpis { //daysnpis lnnpis longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr gen_id) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr idxm) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1920m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1923m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1923m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1910m1, 1926m12), a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/poprobust_nopop_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

* Interpolated population as a control

foreach v in longnpis { //daysnpis lnnpis longnpis {
	
	local treat post_X_`v' pop_int

	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr gen_id) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr idxm) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1920m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1923m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1923m12), a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1910m1, 1926m12), a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/poprobust_popcontrol_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

* 1910 Population

gen pop1910_temp = pop_int if f_myr==tm(1910m4)
bys gen_id: egen pop1910 = mean(pop1910_temp)
drop pop1910_temp

foreach v in longnpis { //daysnpis lnnpis longnpis {
	
	local treat post_X_`v'

	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr gen_id) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr idxm) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1920m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1923m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1923m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1910m1, 1926m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/poprobust_1910pop_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

* 1910 Population with interpolated population as a control

foreach v in longnpis { //daysnpis lnnpis longnpis {
	
	local treat post_X_`v' pop_int

	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr gen_id) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr idxm) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1920m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1920m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1916m1, 1923m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1913m1, 1923m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)
	eststo: ppmlhdfe pat_wtd_inv `treat' if tin(1910m1, 1926m12), exp(pop1910) a(f_myr idxm regyr) cluster(gen_id)

	estout 	using "$RES/fortable/poprobust_1910popcontrol_`v'.tex", replace cells(b(star fmt(%9.3f)) se(par)) ///
		stats(N, fmt(%9.0g) labels("\(N\)"))  starlevels(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
		style(tex) drop(_cons) eqlabels(none) label
	est 	clear
}

